% Simulate data from an SVAR(0) and characterise the likelihood and 
% posterior for an impulse response under a single shock-sign restriction.

clear variables
close all
clc

rng(23032021); % Set random seed

T = 3; % Length of time series
n = 1000; % Number of gridpoints over which to evaluate likelihood
k = 1; % Period in which first structural shock is positive

% Set true values of structural parameters.
A0 = [1 0.5; 0.2 1.2]; % Implies sig21 < 0

% Compute reduced-form parameters
% Lower-triangular Cholesky factorisation of innovation covariance matrix.
Sigmatr = chol(inv(A0)*inv(A0)','lower'); 
Sigma = Sigmatr*Sigmatr'; % Covariance matrix of VAR innovations.
sig11 = Sigmatr(1,1);
sig21 = Sigmatr(2,1);
sig22 = Sigmatr(2,2);

% Construct grid for theta over which to evaluate likelihood
thetaGrid = linspace(-pi+eps,pi-eps,n)';

%% Compute flat regions of likelihood under different realisations of data.
% Appendix to Giacomini, Kitagawa and Read (2022, JBES) presents three
% expressions for the conditional identified set depending on the
% realisations of the data (given sig21 < 0). Consider first two cases
% here:
% (I)  h(phi,y_k) < 0
% (II) h(phi,y_k) > 0 and sig22/sig21 < C(phi,y_k)
% Simulate from DGP until we obtain observations covering these two cases.

% Storage for cases
theta_cis = zeros(2,2);
theta_lik = zeros(2,1);
cases = zeros(2,1);

while ~all(cases) % Keep drawing until exhaust all cases

    % Draw a time series of structural shocks such that the first 
    % structural shock in the kth period is positive.
    epst = randn([T,2]);   
    flag = 0;
    while flag == 0
        epst(k,1) = randn;
        flag = (epst(k,1) >= 0);
    end

    % Use shocks and structural parameters to generate data.
    yt = (A0\epst')';

    % Given simulated data, check which case we're in, compute conditional
    % identified sets for theta, and value of likelihood in flat region.
    h = sig21*yt(k,1)-sig11*yt(k,2);
    C = sig22*yt(k,1)/h;

    if h < 0 % Case (I)

        cases(1) = 1;
        theta_cis(1,1) = atan(max(sig22/sig21,C));
        theta_cis(1,2) = pi + atan(min(sig22/sig21,C));
        theta_lik(1) = exp(log(2) - T*log(2*pi) - (T/2)*log(det(Sigma)) ...
                    - 0.5*trace((yt/Sigma)*yt'));

    elseif h > 0 && sig22/sig21 < C % case (II)

        cases(2) = 1;
        theta_cis(2,1) = atan(sig22/sig21);
        theta_cis(2,2) = atan(C);
        theta_lik(2) = exp(log(2) - T*log(2*pi) - (T/2)*log(det(Sigma)) ...
                - 0.5*trace((yt/Sigma)*yt'));    

    end

end

% Plot likelihoods.
ll1 = (((thetaGrid >= theta_cis(1,1)) & (thetaGrid <= theta_cis(1,2)))*theta_lik(1));
ll2 = (((thetaGrid >= theta_cis(2,1)) & (thetaGrid <= theta_cis(2,2)))*theta_lik(2));
figure;
h1 = plot(thetaGrid,ll1,'color',[0 0 204]./255,'LineWidth',2);
hold on;
h2 = plot(thetaGrid,ll2,'color',[255 0 127]./255,'LineWidth',2);
xlim([-pi pi]);
xticks([-pi -pi/2 0 pi/2 pi]);
xticklabels({'-\pi','-\pi/2','0','\pi/2','\pi'});
xlabel('\theta');
title('Likelihood');
lgnd = legend([h1 h2],...
    {'(I)','(II)'},'interpreter','LaTeX','Location','NorthWest');
legend boxoff;
lgnd.FontSize = 14;
Ax = gca;
Ax.FontSize = 14;
box on; 
cd('Figures');
print('ShockSignLikelihood','-depsc');
print('ShockSignLikelihood','-dpng');
cd ..

%% Approximate posterior for impulse response under uniform prior over theta.
draws = 1e6; % Draws from posterior
etaDraws = zeros(draws,2);

for jj = 1:2

    % Draw theta from its posterior, which is uniform over conditional
    % identified set.
    thetaDraws = theta_cis(jj,1) + (theta_cis(jj,2) - theta_cis(jj,1))*rand([draws,1]);
    % Transform draws into impulse response.
    etaDraws(:,jj) = sig11*cos(thetaDraws);

end

% Plot histogram of posterior for impulse response.
figure;
h1 = histogram(etaDraws(:,1),'Normalization','pdf',...
    'FaceColor',[0 0 204]./255,'EdgeColor',[0 0 204]./255,'EdgeAlpha',0);
hold on;
h2 = histogram(etaDraws(:,2),'Normalization','pdf',...
    'FaceColor',[255 0 127]./255,'EdgeColor',[255 0 127]./255,'EdgeAlpha',0);
xlim([-sig11 sig11]);
xticks([-sig11 0 sig11]);
xticklabels({'-\sigma_{11}','0','\sigma_{11}'});
xlabel('\eta');
title('Posterior Density of $$\eta = \sigma_{11}\cos\theta$$','interpreter','latex');
Ax = gca;
Ax.FontSize = 14;
box on;
cd('Figures');
print('ShockSignPosterior','-depsc');
print('ShockSignPosterior','-dpng');
cd ..

cd('Results');
save('ShockSignIllustration_results.mat');
cd ..
